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We introduce a variant of the Banded Random Matrix ensemble and show, using detailed nu- 
merical analysis and theoretical arguments, that the phonon heat current in disordered quasi-one- 
dimensional lattices obeys a one-parameter scaling law. The resulting /3-function indicates that an 
anomalous Fourier law is applicable in the diffusive regime, while in the localization regime the heat 
current decays exponentially with the sample size. Our approach opens a new way to investigate the 
effects of Anderson localization in heat conduction, based on the powerful ideas of scaling theory. 
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Introduction - Anderson localization, i.e. the complete 
halt of propagation in disordered media due to wave in- 
terference eflFects, is an interdisciplinary field of research 
that addresses systems as diverse as classical, quantum 
and atomic-matter waves. This phenomenon was pre- 
dicted fifty years ago in the framework of quantum (elec- 
tronic) waves by Anderson [1] and its existence has been 
confirmed in recent years by experiments with classical 
[2-10] and matter waves [11, 12]. 

Recently, localization phenomena due to randomness 
have attracted considerable interest in the context of heat 
conduction by phonons [13-15]. A central issue of these 
investigations is the determination of the dependence of 
the heat current J on the system size TV. It has been 
commonly believed that disorder scatters normal modes 
and induces a diffusive energy transport that leads to a 
normal heat conduction described by Fourier's law which 
states that J ~ AT"^. However, many recent studies [14- 
21] suggest that in low dimensional disordered harmonic 
chains this may not always be true. Instead one finds 
that J ~ N~°', where a is usually different from one. 
Although this conclusion is generally accepted for one- 
dimensional systems, where theoretical methods of inves- 
tigation are available, the validity (or not) of Fourier law 
in higher dimensions is totally unclear since the majority 
of the available results are based on numerical simula- 
tions which are limited to small systems sizes. 

In fact, recent experiments on heat conduction in nan- 
otubes and graphene flakes have reported observations 
which indicate such anomalous behaviour with the sys- 
tem size [22-24]. Therefore, not only is it a fimdamental 
demand for the development of statistical physics to un- 
derstand normal and anomalous heat conduction in low 
dimensional systems, but it is also of great interest from 
the technological point of view, since the achievement of 
modern nano-fabrication technology allows one to access 
and utilize such structures with sizes in the range of a 
few nanometers up to few hundred nanometers. 

In this Letter, we approach thermal transport in the 
presence of disorder from a different perspective, namely 
we develop a scaling theory for quasi-one-dimensional 



(ID) random lattices described by a modified Banded 
Random Matrix Ensemble (BRM) . Random matrix mod- 
els played a major role in understanding various proper- 
ties of disordered quantum systems, including the struc- 
ture and statistical properties of their eigenstates [25, 27] 
and eigenvalues [28], the conductance [29], delay times 
[30], etc. Here we introduce a BRM ensemble with band- 
width 26-1-1 that describes an array of coupled oscillators 
with long range (6*'^-neighbor) random couplings, in the 
presence of on-site random pinning which is coupled at 
the left and right edges to a pair of Langevin heat reser- 
voirs. We find that the averaged (rescaled) steady-state 
heat current Jn i^oo) of the phononic excitations for an 
array of size N obeys a one-parameter scaling, i.e. 
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where ^ is a universal function of Jn i^oo) alone, and 
takes the following asymptotic forms 
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with ly ss 0.25. The asymptotic (i.e. TV — > oo) partic- 
ipation number measures the degree of localization 
of the normal modes which dominate the transport. For 
any finite sample of size N the number of these modes, 
/, scales as / ~ N~'^ with 7 « 0.1. The scaling exponent 
of the (actual) heat current Jjv = TV"''' • Jjv ~ TV~" is 
found to be a = 1^+7 ~ 0.35, indicating a violation of the 
Fourier law. Eqs.(l-2) are confirmed in the following via 
detailed numerical simulations, supported by theoretical 
arguments. 

Banded Harmonic Chain Model - We consider a ther- 
mally isolated quasi-one-dimensional harmonic oscillator 

chain with fe'^'^-nearest neighbors coupling. It consists of 
TV equal masses {m — 1) described by the Hamiltonian 
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The corresponding equations of motion are = 
dH/dpn, Pn = —dH/dqn, where qn,Pn are respectively 
the individual oscillator displacements and momenta. 
The last term in Eq. 3 is the harmonic coupling between 
the n-th mass and its b neighbors on the left and right. 
The random spring constants fc„j are chosen to be sym- 



nj 



and uniformly distributed accord- 



metric {knj - n.jnj 

ing to k„j G [-f + 1, f + 1] if < |n - j| < b, and 
h 



nj — otherwise. is a coupling strength parame- 
ter that has to be smaller than 2 and is henceforth set 
to unity. The second term in the Hamiltonian is an on- 
site "pinning" potential with a spring constant e„, ran- 
dom and uniformly distributed in [— ^ -I- 1, ^ -1- l] . The 
offset in these random distributions ensure a positive- 
definite spectrum of the eigenfrequencies, i.e. bounded 
motion of the oscillators. The boundary conditions used 
are qo = qpf+i 0. 
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FIG. 1: (Color online) Temperature (upper) and heat current 
(lower) profiles, for band widths of 6 = 6 (black) and fe = 30 
(red) in a system of size = 1000. The flatness of the heat 
current profile indicates a non-equilibrium steady state. 

Next, we want to study the non-equilibrium steady 
states (NESS) of this chain driven by a pair of Langevin 



(Ornstein-Uhlenbeck) reservoirs set at temperature Tj, 
and Tji respectively, and coupled to the first (last) Nb 
masses with a constant coupling strength A. In all nu- 
merical examples we will set iVf, = 15 and A = 1. The 
coupling to the bath is described by modifying the equa- 
tion of motion for the momentum p„ = —d'H/dqn + 
Er=L,fl (-^Pn + v/2Ai;Cn) wherc 6^ = {Win < 
Nf, otherwise}, 9^ = {1 if n > iV - iV&; otherwise}, 
and Cn{t) is delta-correlated white noise (Cn(OCn' (*')) = 

d„n'S{t-t'). 

The two thermal quantities, local temperature and 
the heat current, can be expressed in terms of elements 
of the covariance matrix C{t) = {x{t) ^ x{t)), where 
X = (gi, . . . 5Ar,pi, . . .pjv)^ is the state vector. Using 
stochastic Ito calculus [13, 31] for the system of Eq.(3) 
we find 



dC 

lit 



= ZC + CZ^+Y, 



(4) 



with the 2N x 2N matrices 
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Y^and Y = ^T^Y^, (5) 



where Y'^ = ^J2n=i^n^N+n and Pj = ej (X) ej is a 
diagonal rank 1 projector for basis vectors (ej)^ — Snj- 
The banded N x A^-matrix K with bandwidth 26 + 1 
encodes all the interactions within the harmonic lattice, 
as described by Eq.(3) 
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The NESS covariance matrix, C°°, can be obtained by 
setting the left hand side of Eq.(4) to zero, resulting in 
the Sylvester equation ZC°° -I- C°°Z^ = -Y. 

The local temperature is simply given by Tn = (Pn) = 
^n+N n+N- To find the expression for the heat current J 
we use the continuity relation dt {Hn) + Jn^ Jn-i — 0, 
where {%„) is the thermal fluctuation average of 



Q'}1 1 ^ ^ 

" 2 E - 9") iPn + Pn+]) - K-j,n {q-n - qn~j) {Pn + Pn-j) , (6) 

ensuring that at any given cross section of the chain, all connections are included (see supplemental material for 
details). In terms of the covariance matrix this yields the expression 

^ n i-\-b 

^ 2 E E ^i-n+bj-n {p°^^j^j^ — C°°^^ + C°°_^^pf — C°°_f_pf) . (7) 
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In the central section of the chain (A''^ < n < N — Ni,), which is not directly coupled to the bath, the NESS heat 
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current has to be independent of n due to continuity, i.e. 
J — Jn- In our proceeding numerics, the bath tempera- 
tures are fixed [37] at 7^ = 2, Tr = f . Additionally, all 
thermal calculations are averaged over f 0^ realizations of 
disorder. An example of the local profiles is displayed in 
Fig.l - in particular, the flat profile seen in J„ confirms 
that continuity is fulfilled. 

Localization Properties - We consider the isolated case 
(A = 0). Substituting (/^(t) = An^i, exp {iuj,^t) results 
in an eigenvalue problem — cj^yly = KA^. Again, note 
that the choice of the random distributions in the banded 
random matrix K ensures positive definite eigenvalues, 
Ldl > 0. The extent of the modes is often characterized 
by their participation number (PN) 
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In Fig. 2, the PN's are plotted versus the eigenfrequencies 
for different finite lattice sizes and a fixed bandwidth of 
6 = 5. All states are localized (P2 < N), yet two win- 
dows are observed: a window of highly localized states 
for higher frequencies and a window of states with larger 
PN for lower frequencies. 




FIG. 2: (Color online) Participation number (PN) versus fre- 
quency for & = 5. Various finite system sizes are delineated 
by colors. Colored stripes indicate the groups of the most ex- 
tended modes, for which P2 > O.GP™'"'. A normalized count 
of selected states yield a scaling / ~ A'^"^ for the integrated 
density of states (IDoS), as shown in the inset, for three rep- 
resentative bandwidths: b = 5 (black •), 6 = 10 (red ■), and 
6 = 12 (blue ♦). The best fit indicates 7 0.1 ± 0.05 (black 
solid line). 

We define the spectral window Suj = Wmax — ^min, 
whereby the modes with the larger PN are supported by 
the condition P2(w) > O.GP^"''. This allows us to find 
a scaling behavior for the integrated density of states 



(IDoS) of these modes 
/ (Slo) 



p{uj)duj ~ N-^, 
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Our results for various values of bandwidth h are shown in 
the inset of Fig. 2. The solid line indicates the best fit with 
7 — 0.1 ± 0.05. The averaged (over the spectral window 
Suj and over disorder realizations) PN is ^tv — (^2)5^, 
reported in Fig. 3 for various fe- values versus Typ- 
ically, more than 10* eigenvectors were used for the av- 
eraging. We find that S,N^oo{b) shows a convergence to- 
wards a finite value Coo(&)- For moderate 6- values this 
asymptotic PN is reached, while for larger bandwidths it 
can be extrapolated from the quotient of two fifth order 
polynomials fitted to the data (dashed/dotted -dashed 
lines in the figure). Inspired by previous studies on the 
locahzation properties of BRM's [25-27], we speculate 
that the asymptotic PN will scale as ^00 ~ b'^ ■ Our ex- 
pectation is nicely confirmed by the numerical data, re- 
ported in the inset of Fig. 3. The best fit indicates that 
rj sa 3.5, and thus ^00(6) ^ b^'^- 
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FIG. 3: (Color online) Main: The disordered and spectral 
(over Siu) averaged PN ^jv against 1/N, for various 6- values. 
Saturation of ^jv as A — >■ 00 is observed for moderate band- 
widths. For larger values of 6 a fitting of the data to a rational 
function of fifth order polynomials (dashed/dot-dashed lines) 
allows to extract the limit . Inset: Parametric dependence 
of Coo on b. The dashed line is the power relation b^'^. 

Scaling Theory - Equipped with knowledge of the lo- 
calization properties of the normal modes of our system 
in Eq.(3), wc now turn to the study of the steady-state 
heat current of Eq.(7). This is formally expressed as 
J = /g p{uj)T{ijj)du! « J^""""^ p{uj)T{uj)du!, where p(w) is 
the density of states and t(w) is the frequency dependent 
transmittance. Since heat is transported significantly 
only by the modes with the larger PN, we confine the 
integration range within the spectral window Sid. These 
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modes have similar localization - and therefore transport 
- properties, as shown previously. Therefore, we approx- 
imate t{uj) by its average value over this window {t)su- 
The remaining integral is then just the IDoS of Eq.(9). 
We can now use knowledge for the transmittance of har- 
monic chains to deduce a scaling relation for the rescaled 
heat current J = J/N~^ oc (t)^^. 

Specifically, the transport theory of disordered media 
predicts that the average transmittance {t)s^ of a dis- 
ordered sample of length N which is characterized by 
a localization length ^oo follows a one-parameter scaling 
{t)suj = /t(A), where the one parameter is A = ^oo /N. 
It is natural to then speculate that the same scaling re- 
lation will apply for the rescaled heat current J. In the 
main panel of Fig. 4 we show our numerics of J plotted 
against A, for a number of different bandwidths and sys- 
tem sizes (6 G [2,35],A^ e [10^,10^]). In this approach, 
wc have used ^oo as a scaling parameter, which allows 
us to collapse all data associated with various 5- values to 
one scaling curve. By visual inspection [38], we find that 
the scaling parameter ^oo ~ ^3-5±o.2^ which confirms the 
previous independent scaling analysis from the participa- 
tion numbers (see Fig. 3 and its discussion). The obvious 
data collapse confirms the conjecture that Jn {£,00) is a 
function of A only i.e. Jn{(,oo) = //(A). 
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FIG. 4: (Color online) Mam: Rescaled heat current J vs. A = 
^oo/N where ~ b'^'^ as was found from Fig. 3. The dashed 
line is a fit to the analytical curve of Eq.(lO). Inset: The 
resulting /3-function of Eq.(l). The dashed lines correspond 
to the asymptotic limits of /? = 1.28 4- 0.94 ■ In J for TV 00 
and /3 = for Af -> 0. 

Next we want to determine the analytical form of the 
scaling function /j(A). We have found that in the limit 
of the localized regime (A ^ 1) this dependence has the 
form J ^ e"*^"/^, in agreement with previous theoreti- 
cal results for pinned harmonic chains with only nearest 
neighbor coupling and mass disorder [16]. In the other 
limit of A 1, the heat transport is diffusive. Assum- 



ing validity of the Fourier law, we may expect a scaling 
of the type J ^ 1/N^~'' - however, recent investigations 
[20, 34] found an anomalous behavior of the heat current, 
which results in the scaling J ^ 1/N°'~^. We therefore 
speculate that in the diffusive domain of A 1, the 
rescaled steady state heat current will follow the relation 
J(A) ~ ciA"^. A possible interpolating law valid in all 
regimes (including the crossover region) is 



J(A) = (c2 + ciA'')exp(-co/A). 



(10) 



Comparison with numerical data in the two limits (A ^ 
1, A < 1) yields v « 0.25, Ci « 0.4 and cq w 0.06. 
Adjusting the last parameter C2 to fit numerics in the 
intermediate region yields C2 ~ 0.012. The resulting an- 
alytical formula nicely fits the numerical results in all re- 
gions, and therefore provides a compact summary of our 
empirical data (see dashed line in Fig. 4). We stress that 
the limiting value of J(A ^ 1) oc l/N^'^'^ leads to an 
anomalous heat exponent a « 0.35 ± 0.05. It should be 
noted that this value is less than what has typically been 
seen in other - nonlinear and disordered - chain models, 
which show values from « 0.4 - 0.7 [18, 20, 21, 34-36]. 

Eq.(lO) can be rewritten in the form of Eq.(l). This is 
the main result of the present Letter, as it allows postu- 
lating the existence of a /^-function for the Jn of generic 
quasi-l-D disordered systems. The resulting /3-function is 
plotted in the inset of Fig. 4. Its asymptotes are seen to 
follow 13 = 1.28 + 0.94 • In J for TV 00, and /S = -v for 
iV ^ 0. 

Conclusions ~ We presented a one-parameter scaling 
theory for the steady-state heat current of quasi-one di- 
mensional disordered harmonic systems with substrate 
pinnings described by a variant Banded Random Matrix 
ensemble. Via numerical analysis and theoretical con- 
siderations, we have established Eq.(l), which allows us 
to conclude that changing disorder strength (or coupling 
range) and system size in the way described by Eq.(lO), 
would not change the renormalized (average) heat cur- 
rent. The one-parameter scaling theory presented here is 
a powerful approach in the quest of understanding ther- 
mal transport, and validity of the Fourier law, in disor- 
dered media. Of further interest will be to investigate 
higher moments of the NESS heat current and also to 
establish a scaling theory for the thermal profile, as a 
function of the scaling parameter A. Although the fo- 
cus of this Letter was on harmonic quasi- ID disordered 
phononic transport, our approach can be used to study 
high-dimensional pinned harmonic systems, and to bet- 
ter understand the effects of phonon-phonon interactions 
[27] in thermal transport. 
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Supplemental Material 



Expression for the Heat Current in Long-Range Connected Oscillator Chains 

The usual way to find the correct expression for the heat current in oscillator chains is to start with the continuity 
relation 

^|^!^ + J„-J„_i=0 (1) 

that connects the local energy density change with the net heat current. The brackets (•) indicate the average over 
the bath stochasticity. By examining the expression for the change in local energy 

—Q^ ^ 2 ^ " ^n) ijPn + Pn+]) - kn-j,n (Qn - 9n-i) {Pn + Vn-j) (2) 

one can identify the terms contributing to J„ and J„_i, respectively. In our chain with longer than nearest neighbor 
coupling, however, this procedure misses terms in the heat current, as illustrated in Fig. 1: Let us examine the current 
Jn at the position indicated by the right dashed vertical line. Apart from the couplings from site n to its neighbors 
on the right (thick lines) there are other couplings penetrating the vertical line. Those connections transport heat 
through this section as well and thus have to be accounted for in the expression for J„. From Eq. 2 we can see, that 
every coupling contributes a term — {ki j{qi — qj)(pi + Pj)) to the current. We thus find 

Jn = - l^ij (in ~ 9j)(P^ + Pj)) ^ (3) 

i>n,j<n 
i-j<b 

that can easily be expressed in terms of entries of the covariance matrix as given in the article. From Fig. 1 we see, 
that the thin lines penetrating the right dashed vertical line penetrate the left one as well, therefore they drop out of 
the expression J„ — J„_i and do not occur in the continuity equation. 
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FIG. 1: Illustration of the contributions to the thermal current for a banded chain model with b=3, i.e. each oscillator is 
coupled to the three neighbors on each side. 



